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Abstract. This article is focused on the study of a convection problem in a 2D set-up in 
the presence of the 0(2) symmetry. In the fluid, viscosity depends on the temperature by 
changing its value abruptly in an interval around a temperature of transition. We explore 
the morphology of the plumes for several parameter settings in the viscosity law and perform 
bifurcation studies at several aspect ratios. We report that at a large aspect ratio and high 
Rayleigh numbers, travelling waves, heteroclinic connections and chaotic regimes are found, 
which are greatly influenced by the presence of the symmetry. 

1. Introduction 

This paper addresses the numerical study of convection at infinite Prandtl number in fiuids 
in which viscosity strongly depends on temperature in the presence of 0(2) symmetry. The 
study of convection in fiuids with temperature-dependent viscosity is of interest because of 
its importance in engineering and geophysics. For instance, there are many fiuids in which 
viscosity varies strongly with temperature; in the Earth's upper mantle viscosity decreases with 
temperature, with contrasts of several orders of magnitude. This problem has been addressed 
both in experiments [42l El E [52] and in theory [35l [40l El [45l |46]. In these contexts, 
the dependence of viscosity with temperature is expressed by means of an Arrhenius or an 
exponential law. Theoretical studies have also treated other dependencies such as linear |37l |4T] 
or quadratic ones [T7[ |49] . 

In this article, we focus on the study of a fiuid in which the viscosity changes abruptly in 
a temperature interval around a temperature of transition. This defines a phase change over a 
mushy region, which expresses for instance the melting of minerals or other components. This 
kind of problem has recently been addressed in [48]. Melting and solidification processes are 
important in the evolution of sea ice and lava lakes [54]; in magma chamber dynamics [53] . 
and in the formation of chimneys in mushy layers [TTJ [TJ [24] , which are also formed in metal 
processing in industry (see, for example, [44]). In phase transitions, other fiuid properties in 
addition to viscosity may change abruptly, such as density or thermal diffusivity. However, in 
this study we consider solely the study of effects due to the variability of viscosity. 

In classical convection problems (with constant viscosity), the study of the solutions and 
bifurcations in the presence of symmetries has been the object of much attention ^1 [23l [36l 
[33l [3T| [30] [T4] [4] , while its counterpart in fiuids with viscosity depending on temperature has 
received less attention. The motivation of this paper is the fact that symmetric systems typically 
exhibit more complicated bifurcations than non-symmetric systems and introduce conditions 
and degeneracies in bifurcation analysis [121 HOI [H] [18] . There exist numerous novel dynamical 
phenomena whose existence is fundamentally related to the presence of symmetry, including 
rotating waves [43 , modulated waves [39l [2], slow "phase" drifts along directions of broken 
symmetry and stable heteroclinic cycles [25l [2 [16 . The S0(2) symmetry is present in 
problems described by the Navier-Stokes [211 Bl] the Kuramoto-Sivashinsky [3l [16] equations 
with periodic boundary conditions, since the equations are invariant under translations and the 
boundary conditions do not break this invariance. Additionally, if the refiection symmetry exists, 
the full symmetry group is the 0(2) group. 

In this context, this paper addresses the convection of a 2D fiuid layer with temperature- 
dependent viscosity and periodic boundary conditions possessing the 0(2) symmetry. Among 
other solutions, we show the presence of travelling waves or limit cycles near heteroclinic connec- 
tions after a Hopf bifurcation, as previously reported in diverse contexts in the presence of this 
symmetry [2] [50] [16] , and here these waves are reported for convection with variable viscosity. 
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Figure 1. Problem set-up. 



Our 2D physical set-up is idealized with respect to realistic geophysical flows occurring in the 
Earth's interior as these are 3D flows moving in spherical shells [71[8]. Under these conditions, 
the symmetry present in the problem is formed by all the orientation preserving rigid motions 
of that fix the origin,which is the S0(3) group [lOl [191 [27] . The link between our simplified 
problem and these realistic set-ups is that the 0(2) symmetry is isomorphic to the rotations 
along the azimuthal coordinate, which form a closed subgroup of SO (3). In this way, the waves 
observed in our setting may be a signature of the presence of travelling azimuthal waves. 

The article is organized as follows: In Section [2| we formulate the problem, providing the 
description of the physical set-up, the basic equations and boundary conditions. In Section 3 
we present the viscosity law under consideration and discuss several limits in which previously 
studied dependencies are recovered. Section 4 summarizes the numerical methods used to sketch 
an outlook of the solutions displayed by the system. Section 5 discusses the solutions obtained 
for a broad parameter set. Finally Section [6] presents the conclusions. 

2. Formulation of the problem 

As shown in Fig [l] we consider a fluid layer, placed in a 2D container of length L [x coordinate) 
and depth d {z coordinate). The bottom plate is at temperature Tq and the upper plate is at 
Ti, where Ti = Tq — AT and AT is the vertical temperature difference, which is positive, ie, 
Ti < To. 

The magnitudes involved in the equations governing the system are the velocity field u = 
{ux^ Uz)^ the temperature T, and the pressure P. The spatial coordinates are x and z and the time 
is denoted by t. Equations are simplified by taking into account the Boussinesq approximation, 
where the density p is considered as constant everywhere except in the external forcing term, 
where a dependence on temperature is assumed, as follows p = po{l — a{T — Ti)) . Here po is the 
mean density at temperature Ti and a the thermal expansion coefficient. 

The equations are expressed with magnitudes in dimensionless form after rescaling as follows: 
{x',z') = {x,z)/d, t' = nt/d^, m' = dM/n, P' = d^P/{poKUo) , 0' = (T-Ti)/(AT). Here, Hi is the 
thermal diffusivity and is the maximum viscosity of the fluid, which is viscosity at temperature 
Ti. After rescaling the domain, Qi = [0, T) x [O^d] is transformed into Q2 = [O^T) x [0, 1] where 
r = L/d is the aspect ratio. The system evolves according to the momentum and the mass 
balance equations, as well a to the energy conservation principle. The non-dimensional equations 
are (after dropping the primes in the fields): 

(1) V-u = 0, 

^fan ^ n . \7^^^ = T^ftPn - \7 P 4- Hiv f ^^^^ ' 

Pr 
dtO- 



-(5^U + U- VU) : 

u • V(9 = Ai9. 



-(Vu+(Vu)^) 



(2) 
(3) 

Here, 63 represents the unitary vector in the vertical direction, R = d'^ag/ST / {v^^n) is the 
Rayleigh number, g is the gravity acceleration, Pr = vq/h\s the Prandtl number. Typically for 
rocks Pr, is very large, since they present low thermal conductivity (approximately l{)~^im? / s) 
and very large viscosity (of the order lO^^TVs/m^) [15 . Thus, for the problem under consider- 
ation, Pr can be considered as infinite and the left-hand side term in ([2| can be made equal 
to zero. The viscosity v{9) is a smooth positive bounded function of ^, which in our set-up 
represents a transition in the fluid, due for instance to the melting of minerals caused by an 
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Figure 2. Representation of the arctangent viscosity law versus the dimension- 
less temperature for different parameters values; a) 6 = 10, a = 0.1 and different 
R values; b) a = 0.1, R = 1300 and different b values; c) 6 = 10, = 1300 and 
different a values. 



abrupt change in viscosity at a certain temperature. This is discussed in detail in the following 
section. 

As regards boundary conditions, we consider that the bottom plate is rigid and that the upper 
surface is non-deformable and free slip. The dimensionless boundary conditions are expressed 
as, 

(4) ^ = 1^ u = 0, on z = and = dzUx = Uz = 0, on z = 1. 

Lateral boundary conditions are periodic. Jointly with equations ([T])-(|3|, these conditions are 
invariant under translations along the x-coordinate, which introduces the symmetry SO (2) into 
the problem. In convection problems with constant viscosity, the reflexion symmetry x 
—X is also present insofar as the fields are conveniently transformed as follows {O^Ux^Uz^p) 
(6>, —Ux^Uz^p). In this case, the 0(2) group expresses the full problem symmetry. The new terms 
introduced by the temperature dependent viscosity, in the current set-up Eq. ^ maintain the 
reflexion symmetry, and the symmetry group is 0(2). 

3. The viscosity law 

We consider that the viscosity depends on temperature, and that it changes more or less 
abruptly at a certain temperature interval centred at a temperature of transition. This is ex- 
pressed with an arctangent law which reads as follows: 

(5) iy{T) = Ai arctan(/3{(T - Ti) - h}) + A2 

The parameter j3 controls how abrupt the transition of the viscosity with temperature is. Very 
high P values imply that the viscosity transition occurs within a very narrow temperature gap, 
while a finite and not too large value j3 assumes that the phase change happens over a mushy 
region of finite thickness [48 . For the results reported in this article, we have fixed /3 = 0.9. As /3 
is fixed, the viscosity transition always occurs in a temperature interval with constant thickness 
~ 0.23. The temperature at which the transition occurs is controlled by h. The constants 
Ai and A2 are adjusted by imposing that at the reference temperature Ti the viscosity law ([5| 
must be vq. On the other hand, in the limit T >> Ti, for instance T — Ti = 2500, the viscosity 
is fixed to a fraction a of the viscosity vq. These conditions supply the system: 

= Ai arctan(— /36) + A2 
uoa = Ai arctan(/3{2500 - b}) + A2 

which has the solution: 

^ ^ ^0(1 -Q^) 

^ arctan(-/36) - arctan(/3(2500 - b)) ' 

^2 = ^0 — ^1 arctan(— 6/3). 
In dimensionless form, the viscosity law becomes: 

(6) ^ = Ci arctan(/3(i?l9/i - b)) + C2 

where Ci = Ai/uq and C2 = A2/iyo- In this expression, R is the Rayleigh number, is the 
dimensionless temperature, which takes values between at the upper surface and 1 at the 
bottom. The parameter /i, defined as /i = v^n/ {d^ag)^ is in this study fixed to /i = 0.0146. 
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Figure 3. The law of the viscosity dependent on the temperature used in 
p8] with viscosity contrast of factor 10 against the arctangent law (|6| with 
parameters b = 1, 5, 10, 30, R = 2500 and a = 0.1 



The parameter a is related to the inverse of the maximum viscosity contrast on the fluid layer, 
although the viscosity u^a may not correspond to any element of the fluid layer. For instance 
Figure |2^) shows the viscosity variation with temperature for different Rayleigh numbers at 
a = 0.1 and b = 10. It is observed that at low R = 600, the viscosity is almost uniform 
in the fluid layer, and it is only beyond R = 1000 that the sharp change in the viscosity is 
perceived. Figure [2]3) shows the effect of varying 6 at = 1300 and a = 0.1. If b is as small as 
1, the transition occurs close to ^ = 1 and most of the layer has low viscosity, while if b is very 
large at this R number most of the fluid has constant viscosity uq. It is interesting to relate the 
viscosity law as represented in these figures with the linear stability analysis of a fluid layer with 
constant viscosity uq^ as presented in Figure |4] In this figure, one may observe that the critical 
R number is approximately Rc ~ 1100. On the other hand, in Figure ^p) one may observe 
that if b is large, the viscosity near the critical Rayleigh number is almost constant across the 
fluid layer. In this case, the phase transition is noticed in the fluid at large R numbers, well 
above R = 1300, in a convection state in which vigorous plumes are already formed, as may 
be deduced from Figure [2|i). Figure [4^) confirms that at this limit the instability threshold of 
the conductive state remains unchanged with respect to that obtained with constant viscosity. 
On the other hand, if b is small, changes in the fluid viscosity are noticed at low R numbers 
-below the critical threshold of a fluid with constant viscosity- and in this case the instability 
threshold of the conductive state is affected by the phase transition. This is illustrated, for 
instance, in Figures |2|l) and|4]3). For b = 10 and a = 0.1, the changes in the viscosity across the 
fluid layer are noticed from R = 800 onwards, which is below the instability threshold obtained 
for constant viscosity. In this case, the instability thresholds for the conductive solution are as 
those displayed in Fig|5) and thus the phase transition is perceived from the beginning by weakly 
convective states. 

We now discuss the relation between the arctangent law and an Arrhenius type law frequently 
used in the literature to model mantle convection problems. This viscosity law is expressed 
according to 011 US] as: 



where is the activation energy, R is the universal gas constant, AO is the temperature drop 
across the fluid layer and ti is the surface temperature divided by the temperature drop across 
the layer. Figure 3 represents the viscosity ill versus the dimensionless temperature for = 
0.25328 and ti = 0.1 as considered by [28]. Additionally, several arctangent laws with different 
b values are displayed. In this representation, one may observe the great similitude between the 
Arrhenius law and the arctangent law for 6 = 1. At larger b values, the decaying rate between 
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Figure 4. Critical instability curves i?(m, F) for a fluid layer a) with constant 
viscosity; b) with temperature dependent viscosity ja = 0.0146 a = 0.1 and 
6 = 30. 




Figure 5. Critical instability curves i?(m, F) for a fluid layer with temperature 
dependent viscosity /i = 0.0146 a = 0.1, 0.01 and b = 10 



viscosities is still similar to an Arrhenius law; however, there are temperature intervals exist 
with approximately constant viscosities and uoa. 

4. Numerical methods 

Analysis of the solutions to the problem described by equations ([l])-([3| and boundary condi- 
tions Q is assisted by time dependent numerical simulations and bifurcation techniques such 
as branch continuation. As highlighted by [131 [38], the combination of both techniques provide 
a thorough insight into the solutions observed in the system. A full discussion on the spectral 
numerical schemes used is given in [13 . For completeness, we now summarize the essential 
elements of the numerical approach. 

4.1. Stationary solutions and their stability. The simplest stationary solution to the prob- 
lem described by equations ([l])-([3| with boundary conditions Q is the conductive solution which 
satisfies Uc = and Oc = —z + 1. This solution is stable only for a range of vertical temperature 
gradients which are represented by small enough Rayleigh numbers. Beyond the critical thresh- 
old Rc, a convective motion settles in and new structures are observed which may be either time 
dependent or stationary. In the latter case, the stationary equations, obtained by cancelling the 
time derivatives in the system ([l])-(|3| are satisfied by the bifurcating solutions. As in the con- 
ductive solution, the new solutions depend on the external physical parameters, and new critical 
thresholds exist at which stability is lost, thereby giving rise to new bifurcated structures. These 
solutions are numerically obtained by using an iterative Newton-Raphson method. This method 
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starts with an approximate solution at step s = 0, to which is added a smah correction in tilda: 

(8) (u^ + £i,r +^,P^+P). 

These expressions are introduced into the system ([T])-([3|, and after cancehing the nonhnear terms 
in tilda, the following equations are obtained: 

(9) o=v•a + v•u^ 

^0 

(10) +Ll3(^^)^. + Pl4(^^<,<)^~], 

= -d,P- d,P' + -[P2l(^^ <, Ul) + L22{0nUa: 
^0 

(11) ^L23mu, + (L24(r,<,0 +i?)^], 

(12) =u ■ Vr + • + • - - AO'. 

Here, Lij {i = 1,2, j = 1,2,3,4) are linear operators with non-constant coefficients, which are 
defined as follows: 

(13) Lii{0, Uj,, Uz) =2deiy{0)dx0dxUx + iy{0)Aux + deiy{0)dzO{dxUz + d^u^), 

(14) Li2{0) =2d0u{O)dJd, + u{0)A + d0u{O)d,ed,, 

(15) Li3(^) =S,^(^)S,^S„ 

Li4(6>, Uz) =2dey{9)dxUxdx + 2dlQy{9)dx0dxUx + deiy{0)Aux 

(16) + (5,^/, + S,ii,)(S,^(^)S, + S,2,^(^)S,^), 

(17) L2i{0, u,, u,) =2d0iy{O)d,Od,u, + u{0)Au, + S^^(^)S,^(S,ii, + S,ii,), 

(18) L22{0) =d0u{O)dJd,, 

(19) L23(^, ^x., ^i.) =25^^(^)S,^a, + ^(^) A + d0iy{O)dJd,, 

L24(6>, ^la;, ^^) =2dev{0)dzUzdz + 2dee^{0)dz0dzU^ + dev{0)Au^ 

(20) + (S,^i, + S,ii,)(S^i/(^)S, + S^^i/(^)S,^). 

The unknown fields u, P, ^ are found by solving the linear system with the boundary conditions: 

(21) ^ = 0, u = 0, on 2: = and = dzUx = '^2 = 0, on z = 1. 
Then the new approximate solution s + 1 is set to 

^.+1 = + a, = r + ^, p"+^ = p' ^P. 

The whole procedure is repeated for 5 + 1 until a convergence criterion is fulfilled. In particular, 
we consider that the P norm of the computed perturbation should be less than 10~^. 

The study of the stability of the stationary solutions under consideration is addressed by 
means of a linear stability analysis. Now perturbations are added to a general stationary solution, 
labelled with superindex b: 

(22) u{x,z,t) = u^(x,z) + u(x,z)e^*, 

(23) 0{x,z,t) = 0\x,z) ^0{x,z)e^\ 

(24) P{x,z,t) = P\x,z)^P{x,z)e^K 

The sign in the real part of the eigenvalue A determines the stability of the solution: if it is 
negative, the perturbation decays and the stationary solution is stable, while if it is positive the 
perturbation grows over time and the conductive solution is unstable. The linearized equations 
are: 

(25) =v • a 

(26) = - S,P + - [Pi2(^')^. + L,s{0')u, + Pl4(^^ ul ul)0] 

^0 

(27) = - S,P + - [P22(^')^. + L23{0')u, + (P24(^', ul u^) + R)0] 

^0 

(28) =a • V(9^ + • + • V(9^ - AO ^ XO, 
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where the operators Lij are the same as those defined in Eqs. (13)-(20). Equations (25)-(28) 



jointly with its boundary conditions (identical to (21)) define a generalized eigenvalue problem. 
The unknown fields Y of the stationary ([9l)-(|12[) and eigenvalue problems (25)-(28) are ap- 



proached by means of a spectral method according to the expansion: 

\L/2] M-1 \L/2] M-1 

(29) r(x, ^) = E E ^TmTmiz) cos((/ - l)x) + ^ E ^L^m(^) sin((/ - l)x). 

1=1 m=0 1=2 m=0 

In this notation, [•] represents the nearest integer towards infinity. Here L is an odd number 
as justified in [13^. 4 x L x M unknown coefficients exist which are determined by a colloca- 
tion method in which equations and boundary conditions are imposed at the collocation points 

27r 

Uniform grid: Xj = {j — 1)—— , j = 1, . . . , L; 

Gauss-Lobatto: Zi = cos (^(^J^ — ~ ^) ^) ' ^ = 1, • • • , M; 

according to the rules detailed in |13j. Expansion orders L and M are taken to ensure accuracy 
on the results: details on their values are provided in the Results section. 

4.2. Time dependent schemes. Together with boundary conditions Q, the governing equa- 
tions ([l])-([3| define a time-dependent problem for which we propose a temporal scheme based 
on a spectral spatial discretization analogous to that proposed in the previous section. As be- 
fore, expansion orders L and M are such that they ensure accuracy on the results and details 
on their values are given in the following section. To integrate in time, we use a third order 
multistep scheme. In particular, we use a backward differentiation formula (BDF), adapted for 
use with a variable time step. The variable time step scheme controls the step size according to 
an estimated error E for the fields. The error estimation E is based on the difference between 
the solution obtained with a third and a second order scheme. The result of an integration at 
time n + 1 is accepted if E is below a certain tolerance. Details on the step adjustment are found 
in [13]. 

BDFs are a particular case of multistep formulas which are implicit, thus the BDF scheme 
implies solving at each time step the problem (see [26^): 

(30) = V • u^+^ 

(31) = i^r+^es - VP^+^ + div f ^^^^^(Vu"+i + (Vu^+^)^) 

V ^0 

(32) = -u^+^ • Vi9^+^ + Ai9^+\ 
where dtO^^^ is replaced by a backward differentiation formula. 



In p!3], it has been proved that instead of solving the fully implicit scheme (30)- (32), a semi 



implicit scheme can produce results with a similar accuracy and fewer CPU time requirements. 



The semi-implicit scheme approaches the nonlinear terms in Eqs. (30)-(32) by assuming that the 
solution at time n + 1 is a small perturbation Z of the solution at time n; thus, z"^+^ = z"^ + Z. 
Once linear equations for Z are derived, the equations are rewritten by replacing Z = z^+^ — z^. 
The solution is obtained at each step by solving the resulting linear equation for variables in 
time n -h 1. 

5. Results 

5.1. Exploration of stationary solutions in the parameter space. In this section we 
explore how stationary solutions obtained at a low aspect ratio F = 3.4 for the system ([l])-([3| 
depend on the parameters a and h of the viscosity law (|6|. We examine the shape and structure 
of the plumes in a range of Rayleigh numbers from R = 2500 io R = 3500. 

We first consider that the parameter h is large: for instance, as large as 30. In this case. 
Figure |2]3) confirms that at the instability threshold the viscosity across the fluid layer is almost 
constant and equal to z/q, no matter what the value of a may be. Thus, the viscosity transition 
becomes evident in the fluid once convection has settled in at R numbers well above the instability 
threshold. Figure|6|i) shows the plume pattern observed dX R = 2500 for a = 0.1; although values 
a = 0.01 and a = 0.001 are not displayed, they provide a very similar output. The plume is 
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Figure 6. Plumes obtained for the viscosity parameter b = 30. a) R = 2500 
and a = 0.1; h) R = 3500 and a = 0.1; c) R = 3500 and a = 0.001. 




Figure 7. Plumes obtained for the viscosity parameter b = 10. a) R = 2500 
and a = 0.1; h) R = 3500 and a = 0.1; c) R = 2500 and a = 0.001. 



spout-shaped, with the tail of the plume nearly as large as the head. In the pattern, the two 
black contour lines mark temperatures between which the viscosity decays most rapidly. These 
correspond to the transition region in which the gradient of the viscosity law (|6| is large. Thus 
one of the contours, the coldest one, fits the temperature Oi at which the viscosity has decayed 
by 5% from the maximum, i.e., v = 0.95 i^o, while the second addresses 62 = Oi -\- AO with 
temperature increment = 0.23. The maximum viscosity decay rate always takes place at a 
constant temperature increment, since the decaying rate of the law (61) /3, is the same through 
out all this study. At larger Rayleigh numbers, R = 3500, Figure [63) shows that the head of 
the plume becomes more prominent. A comparison between Figure |6p) and Figure [g]^) indicates 
that the large viscosity contrast favours the formation of a balloon-shaped plume, with a thinner 
tail and more prominent and rounded head. As regards the velocity fields, none of these patterns 
develop a stagnant lid at the surface for any viscosity contrast a, even though the upper part 
corresponds to the region with maximum viscosity. 

We now consider that the parameter b is small. As explained in Section 3, in this case the 
viscosity transition occurs at low R numbers, below the instability threshold of the fiuid with 
constant viscosity uq. As low viscosity also implies diminishing the critical R number, the overall 
effect is that for small b the instability threshold is below that with constant viscosity and 
the phase transition is perceived by weakly convective states. Figure [7^) shows the structure 
of the plume obtained for b = 10 and a = 0.1 at = 2500. The head tends to be spread in 
a wide area and the viscosity transition occurs at cold fiuid zones away from the main plume. 
This pattern is rather similar to those obtained with 6 = 5 or 6 = 1, except that for smaller b 
values the tail of the plume tends to be thinner. Increasing the R number makes the tail of the 
plume thinner and spreads the head of the plume in the upper part, as refiected in Fig. [7]3). 
On the other hand, high R numbers shift the viscosity transition towards colder temperature 
contours. As expected from the viscosity law ([6|, there is no R number at which the whole fiuid 
layer is "melted", since this law always imposes that a transition occurs across the fiuid layer. 
Fig. [t]^) reports the effect of diminishing the viscosity contrast a to a = 0.001 at R = 2500. A 
mushroom-shaped plume with a thin tail and prominent head is observed. As before, none of 
these solutions develop a stagnant lid at the surface for any of the examined viscosity contrasts 
a. 

Intermediate values such as 6 = 17 interpolate these extreme patterns. Fig. [S^i) shows the 
evolution from Fig. [T^) to Fig. |6|i) in which the black contour lines indicating the position of 
maximum viscosity decay, converge towards the ascending plume boundary, thus highlighting its 
shape. The head of the plume shrinks and the tail strengthens. Diminishing a to the contrast 
0.001 transforms the structure into a balloon-shaped plume (Fig. [sfD)), while an increase in the 
R number spreads the head of the plume in the upper fiuid towards a mushroom-shaped plume. 

The results reported in this section are obtained with expansions (L = 37 x M = 44) except 
that in Fig. 7 c), which corresponds to (L = 47 x M = 42). Similarly to what is reported in 
[13]. The validity of these expansions is decided by ensuring that it provides accuracy in the 
eigenvalue along the neutral direction due to the SO (2) symmetry, which is always 0. 
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Figure 8. Plumes obtained for the viscosity parameter b = 17. a) R = 2500 
and a = 0.1; h) R = 2500 and a = 0.001; c) R = 3000 and a = 0.001. 




Figure 9. Bifurcation diagram as a function of the aspect ratio at = 1300 
for a fluid with viscosity dependent on temperature (6 = 10, a = 0.1). 



5.2. Bifurcation diagrams and time dependent solutions. Solutions to the system ([T])-(|3| 
experience bifurcations depending on the aspect ratio and on the Rayleigh number. We now 
describe how these solutions vary along the dotted lines enhanced in Figure [5] for parameters 
/i = 0.0146 and b = 10. We consider for a the choices 0.1 and 0.01. 

Figure [9] shows the branch bifurcation diagram as a function of the aspect ratio for Rayleigh 
number R = 1300 and a = 0.1. Branches are obtained by representing along the vertical axis 
the sum of the absolute value of two relevant coefficients in the expansion of the temperature 
fleld, bii and 6?2- Solid lines stand for stable branches, while dashed lines are the unstable ones. 
At a low aspect ratio, the stable branch is that with wave number m = 1, and at a higher aspect 
ratio the stable solutions increase their wave number to m = 2 and m = 3. The unstable branch 
ending up with a saddle-node bifurcation and connecting the m = 1 with the m = 2 branch 
corresponds to a mixed mode. 

Stationary stable and unstable solutions, obtained at the positions indicated by arrows, are 
pictured. In the patterns, the two black contour lines mark temperatures between which the 
viscosity decays most rapidly. No stagnant lid appears at the surface for any of the aspect ratios 
considered. The expansion orders required by this flgure to ensure accuracy are not the same 
along all branches. A rule of thumb is that high modes obtained at larger aspect ratios require 
higher expansions. Thus while for mode m = 1 expansions (L = 37 x M = 44) are sufficient, for 
m = 2 and m = 3 at larger aspect ratios expansions are increased up to (L = 61 x M = 44). 

Bifurcations are further analyzed at three different aspect ratios as a function of the Rayleigh 
number. Figure 10 represents the branching obtained at F = 3.4 for a = 0.1. The pictured 
plumes, which are computed for a rather low Rayleigh number, R = 1500, are spout-shaped, 
with the tail of the plume nearly as large as the head. As already reported in the previous section 
for increasing R numbers, plumes become balloon-shaped and beyond that mushroom-shaped. 
No stagnant lid is observed at any R number. The black contour lines highlight the temperatures 
between which the viscosity decays more rapidly. Several branches are distinguished. The branch 
related to mode m = 1 arises at the lowest R number and is stable in the whole range displayed. 
Mode m = 2 emerges at ~ 860 from the unstable conductive solution through an unstable 



10 



J. CURBELO AND A.M. MANCHO 




Figure 10. Bifurcation diagram as a function of the Rayleigh number for a 
fluid with viscosity dependent on temperature {b = 10, a = 0.1) at F = 3.4. 



branch, which becomes stable through a pitchfork bifurcation at ~ 890. Results at this aspect 
ratio are obtained with expansions (L = 37 x M = 44). 

This simple diagram with simple stationary solutions obtained at a low aspect ratio is in 
contrast to those with more complex solutions obtained at a larger aspect ratio. Figu re pT] 
represents the bifurcations obtained at F = 6.9 as a function of R for a = 0.01. Figure |11^) 
examines the R interval from 800 to 1300. In this range several stationary solutions are portrayed 
both stable and unstable. As in previous diagrams black contour lines mark the temperatures 
at which the viscosity has the largest decaying rate. At R ^ 1290, a Hopf bifurcation occurs at 
the branch of mode m = 3. After the bifurcation, a travelling wave is found, as illustrated in 
the phase portrait represented at = 1300. The solution evolves in time by travelling towards 
the left. This breaks the symmetry x —x. However, the right travelling solution obtained 
by the symmetry transformation also exists, as expected from equivariant bifurcation theory 
[12]. The presence of travelling waves after a Hopf bifurcation has been reported in diverse 
contexts in under the presence of the 0(2) symmetry [H [50l [161 |T2] , here they are reported 
in the context of convection with variable viscosity. At larger R numbers, up to ~ 1320, the 
travelling wave persists, while its frequency increases. A stable fixed point with wavenumber 
m = 3 is found in the range R ^ 1340 — 1380. A cycle limit appears at around R ~ 1400. In this 
regime, the time-dependent solution consists of plumes that weakly oscillate in the horizontal 
direction around their vertical axis of symmetry. Close to R ^ 1416, a stable branch of fixed 
points emerges, which is visualized at ^ 1525. It shows the presence of plumes that are non- 
uniformly distributed along the horizontal coordinate: two close plumes, which are asymmetric 
around their vertical axis, and a third one that maintains its symmetry. None of the described 
solutions develop stagnant lids at the surface. At low R numbers {i.e. Figure 111)) results are 
obtained with expansions (L = 47 x M = 44), while for higher R numbers {i.e. Figure pTjp)) 
results are obtained with expansions (L = 61 x M = 44). 

Figure 12 shows the bifurcation diagram obtained at F = 7.4 as a function of R for a = 0.1. 
The mode m = 3 branch, marked with a solid black line, emerges at R ^ 794. At R ^ 2190 
the branch undergoes a Hopf bifurcation. Beyond this point, solutions embedded in a projection 
over the coefficient space are represented at the R values marked with vertical dotted lines. A 
limit cycle is observed at = 2210 just above the bifurcation point. Its projection over the 
coefficient space displays a point at every time step of the time series. The solution appears to 
reside in the neighbourhood of a heteroclinic connection between two fixed points as it evolves 
into a quasi-stationary regime -near the large density of points- followed by a rapid transition 
to a new quasi-stationary regime. The two fixed points between which the solution oscillates are 
similar to the non-uniformly distributed plumes described in the previous paragraph. A solution 
is found at R = 2300 that has a time-dependence in which the block of plumes shifts irregularly 
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Figure 11. Bifurcation diagrams as a function of the Rayleigh number for a 
fluid with viscosity dependent on temperature {b = 10, a = 0.01) at F = 6.9. a) 
Rayleigh number in the range 800-1300; b) Rayleigh number in the range 1250- 
1500. 



along the horizontal direction, towards both the left and the right. For increasing R numbers, the 
horizontal motion persists, but the oscillation becomes more regular and pattern displacements 
along the x— coordinate are gradually reduced. This is verified through simulations at = 2350 



and at R = 2400. The diagram displayed in Fig. 12 1) shows a gray solid line associated to a 
mode m = 2 stable branch that emerges by means of a saddle node bifurcation jointly with an 
unstable branch. An irregular pattern obtained at = 1800 for the unstable branch is included 
in this diagram. Once again, none of the solutions described at this aspect ratio has a stagnant 
lid at the surface. Results in this figure are obtained with different order expansions. At low 
R number expansions (L = 47 x M = 42) are sufficient while for higher R numbers they are 
increased up to (L = 61 x M = 44) and even to (L = 101 x M = 44). 

6. Conclusions 



This article addresses the study of a convection problem with temperature-dependent viscosity 
in the presence of the 0(2) symmetry. In particular, the considered viscosity law represents a 
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Figure 12. Bifurcation diagrams as a function of the Rayleigh number for a 
fluid with viscosity dependent on temperature (6 = 10, a = 0.1) at F = 7.4. a) 
Rayleigh number in the range 700-1800; b) Rayleigh number in the range 1800- 
2500. 



viscosity transition at a certain temperature interval around a temperature of transition. This 
is a problem of great interest for its many applications in geophysical and industrial flows. 

Our results report the influence on parameters a and h of the viscosity law on the morphology 
of the plumes at a low aspect ratio (F = 3.4). It is shown that if the temperature of transition 
is well above the instability threshold of a fluid with constant viscosity z^o, ^-e, h is large, plumes 
tend to be thicker and show spout-like shapes. Increasing the R number induces their evolution 
towards balloon-shaped plumes, and this effect is more pronounced for high viscosity contrasts 
(small a). At low h values plumes are thinner, and the head of the plume tends to spread in a 
mushroom-like shape in the upper part of the fluid. 
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We explore bifurcations both for a fixed R number as a function of the aspect ratio, and 
bifurcations at three fixed aspect ratios as a function of the R number. No stagnant hd regime 
is observed in any of the physical conditions analyzed. Among the stationary solutions obtained 
along the bifurcation branches, one of the more interesting stable patterns consists of the non- 
uniformly distributed plumes that break symmetry along their vertical axis. 

We also find that, for the higher Rayleigh numbers explored, at a high aspect ratio several rich 
dynamics appear. As already reported in classical convection problems, we find dynamical phe- 
nomena fundamentally related to the presence of symmetry, such as travelling waves, oscillating 
solutions in the neighbourhood of heteroclinic connections and chaotic regimes characterized by 
"phase" drifts along the horizontal direction linked to the SO (2) symmetry. 
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